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Abstract 

We show a numerical scheme to solve the moment equations of the radiative 
transfer, i.e., Ml model which follows the evolution of the energy density, E, and the 
energy flux, F. In our scheme we reconstruct the intensity from E and F so that 
it is consistent with the closure relation, relation, x = (3 + 4/^)/(5 + 2a/4 — 3/^). 
Here the symbols, Xi f = and c, denote the Eddington factor, the reduced 

flux, and the speed of light, respectively. We evaluate the numerical flux across the 
cell surface from the kinetically reconstructed intensity. It is an explicit function of E 
and F in the neighboring cells across the surface considered. We include absorption 
and reemission within a numerical cell in the evaluation of the numerical flux. The 
numerical flux approaches to the diffusion approximation when the numerical cell itself 
is optically thick. Our numerical flux gives a stable solution even when some regions 
computed are very optically thick. We show the advantages of the numerical flux 
with examples. They include flash of beamed photons and irradiated protoplanetary 
disks. 
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1. Introduction 

Radiation plays important roles in many astronomical objects and media. Thus we need 
to include the effects of radiation somehow in a realistic numerical simulation. However it is 
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still an extremely heavy load to solve the full radiative transfer, i.e., taking account of energy 
spectrum and angular distribution as well as spatial distribution and temporal evolution. This 
is simply because the radiative intensity is a function in 6 dimensional phase space. We need 
to reduce the load by using some approximations. 

There exists several types of ideas for reducing the computation cost. Each of them have 
advantages and weak points. 

First we can use a symmetry to reduce the cost for computation. If the spherical sym- 
metry is a good approximation, the intensity is a function in three dimensional phase space 
and the radiative transfer is relatively easy to solve. However this idea can be applied only to 
highly symmetric systems. 

Second we can assume that only a few sources are dominant in the radiation fields. If 
most of light comes from restricted relatively small number of sources, solving radiative transfer 
is relatively easy. We need to solve only the light rays from the sources. The number of light 
rays to be solved are greatly reduced. It is a good approximation when some stars or black 
holes are dominant light sources (see, e.g., Susa 2006; Okamoto, Yoshikawa & Umemura 2012, 
and the references therein). However reflection and reemission from diffuse media cannot be 
taken into account under this approximation. 

Third we can use the moments of radiative intensity such as the energy density by using 
some approximation on the angular distribution such as the flux limited diffusion (FLD, see 
e.g., Levermore 1984). FLD uses only the radiative energy density, which is equivalent to the 
average radiative intensity over the solid angle of photon traveling, to express the radiation 
field. The radiative flux is derived from the gradient of the energy density. This approximation 
is good when the medium is optically thick and radiation field is nearly isotropic. It is not 
so bad even when the medium is transparent. However, this approximation cannot express a 
shadow, a dark region behind an opaque object of finite extent. Spurious radiation erases out 
the shadow in FLD (Gonzalez et al. 2007). Scattering is another weak point of FLD. Scattering 
changes angular distribution of intensity but not the energy density. Thus FLD cannot evaluate 
the effects of scattering since FLD takes account of only energy density. 

Gonzalez et al. (2007) proposed Ml model in which the radiation field is expressed by 
the energy density and flux, i.e., both the 0th and 1st moments of the radiative intensity. 
It is demonstrated that Ml scheme can simulate a shadow by an opaque sphere successfully. 
Ml scheme can take account of scattering. Scattering reduces the energy flux while keeping 
the radiation density constant. When the scattering is isotropic, the effect is correctly taken 
into account in Ml model. However, also Ml model has several weak points and limitations. 
Ml model cannot solve crossing of two beamed lights. They erroneously merge into a beam 
at the crossing point. This limitation comes from the fact that Ml model evaluates higher 
moments of radiation intensity from the 0th and 1st moments. In other words. Ml model has 
the lowest angular resolution and crossing of two beams are beyond the scope. This limitation 
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is compensated by the low computation cost. 

Ml model equations are similar to the hydrodynamical equations in the conservation 
form. Both the equations are hyperbolic and have source terms. Thus we can apply numerical 
methods for integration of the hydrodynamical equations to solve Ml model equations. However 
we have two concerns when solving Ml model equations. First the characteristic speeds are 
complex and not easy to evaluate, though modern schemes for numerical hydrodynamics rely 
on them (see, e.g., Toro 2009). We can avoid the computation of characteristics by using HLLE 
flux but the resultant flux is diffusive and makes a shadow dim as pointed out by Gonzalez et 
al. (2007). 

Second, absorption and emission (the source terms) are dominant when optically thick. 
On the other hand, the source terms due to gravity are minor contributions in numerical 
hydrodynamics. Thus they are simply added after solving wave propagation. This approach 
does not work well in Ml model equations when a cell itself is optically thick. This difficulty is 
known as the diffusion limit behavior and several solutions are proposed in the literature (see. 
e.g.. Audit et al. 2002; Berthon, Charrier, & Dubroca 2007). 

In this paper we propose an idea to construct a numerical flux for Ml model which 
is less diffusive and yet stable in the diffusion limit. First we show a method to evaluate a 
numerical flux of Ml model from radiation intensity kinetically reconstructed from the radiation 
energy density and flux. The reconstructed radiation intensity is consistent with the closure 
relation, i.e., the formula to close the moment equations of the radiative transfer. We evaluate 
the radiative flux and pressure across the boundary between two adjacent computation cells 
by integrating the reconstructed intensity over the solid angle. We use the intensity of the 
upwind side, the numerical flux is subject to causality. Fortunately, the numerical flux is an 
explicit function of the radiation energy density and flux. A similar scheme is constructed for 
gas djTiamics and called "kinetic scheme" (see, e.g., Pullin 1980; Deschpande 1986 and the 
references cieted in Hauck 2011). Thus we use the same terminology in this paper. 

Second we include absorption within a computation cell. The numerical flux is evaluated 
on the cell surface, while the energy density and flux are evaluated at the cell center. They can 
be appreciably different if the cell itself is optically thick. We propose an interpolation formula 
which provides a good approximation both in the optically thin and thick limits. It approaches 
to one obtained from the reconstructed intensity in the optically thin limit, while it does to the 
diffusion approximation in the optically thick limit. We show that this numerical flux gives a 
stable solution even when the computation box contains both optically thin and thick cells. 

This paper is organized as follows. We describe our numerical methods to solve the Ml 
model in §2. We perform some simple examples to show the nature of our numerical scheme 
in §3. In §4, we apply Ml model to irradiated protoplanetary disks. We discuss accuracy and 
stability of our numerical flux in §5. We also discuss affinity of Ml model for massively parallel 
computing in §5. Methods for constructing numerical flux of the second order accuracy in space 
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are given in Appendix. 



2. Ml Model 

2.1. Basic Equations 

First we review Ml model of Gonzalez et al. (2007). We assume that emission is thermal 
and scattering is isotropic. Then the radiative transfer equation for the specific intensity, 
Iu{x, t; n), is expressed as 

(-^ + n-V]h{x, t; n) = Ky^aByp{x, t) - {n^^a + i^u,s) plu{x, t; n) 



c dt 



+ f^u,sP J luix, t; n') dn' , (1) 



where p and c denote the density the speed of light, respectively; hiu,a and k^^s denote the 
absorption and scattering opacities at the photon frequency, z/, respectively. The symbols, n 
and n', denote the angular variable, i.e., the unit vector parallel to the light propagation. The 
symbol, B,^, denotes the Planck function and is a function of the temperature, T. 
We integrate Equation (1) over the solid angle to obtain 

BE ^ BE ■ 

" ■ ElT^ = ^u,aP{^nB, - cE^) , (2) 



dt ^ dxi 



where 



Ey^i{x, t) = J (ei-n) I^{n)dn. (4) 

The symbol, i, specifies a direction in the Cartesian coordinates and does the unit vector 
in the direction. Equation (2) denotes the conservation of radiation energy density, E^, when 
the right hand side vanishes. For simplicity we assumed that the emission and absorption are 
isotropic. 

Similary we obtain 

~~Q^ ^ ^ = - C {Ku,a + Hu,s) P Fu,i , (5) 

where 

P^^ij {x, t) = ^ J {si- n) {sj ■ n) Iy{n) dn , (6) 

by integrating Equation (1) multiplied by n over the whole solid angle. The symbol, j, as well 
as i specifies one in the Cartesian coordinates. 

When deriving Equation (5), we assumed that the scattering is symmetric with respect 
to forward and backward. When the scattering is anisotropic. Equation (5) should be replaced 
by 
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dt 



3 f)p .. 

^ dx- 



+ — (cos 9))] pFy 



(7) 



where (cos 9) denotes the scattering asymmetry parameter. 

In order to solve Equations (2) and (5), we invoke the closure relation, 



\fu\ 



where 



fu,2 
V fu,3 J 
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cE, 



\ J 



(8) 



(9) 



(10) 



5 + 2^4-31/,, 

The above radiative transfer equation is solved in coupled with the hydrodynamical 
equations. The solution of the hydrodynamical equations provides the density, velocity, and 
temperature. Hence the opacity and source functions are evaluated as the functions of them. 
In this paper we consider only the change in the temperature and neglect the changes in the 
density and velocity. This approximation can be justified when we consider the protoplanetary 
disk in thermal equilibrium. The frequency dependent opacity depends little on the density and 
temperature in a certain regime (see, e.g., Henning & Stognienko 1996), although the Rosseland 
mean opacity does depend. 

2.2. Hydrodynamics 

The gas is heated by absorption and cooled by emission. The heating and cooling are 
evaluated by 

Ds 



Us f°° 
pT— = / a,y^a[cE^ - 47rB^(T)]du, 
Ut Jo 



(11) 



where D/ Dt and s denote the Lagrange derivative and the specific entropy, respectively. Then 
the hydrodynamical equations are written in the conservation form as 



dU 



H 



dt 



where 



E 

i=l 



dF 



HA 



dxi 



(12) 
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pV2Vi + P62,i 

pHuVi 



(13) 
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P9i 

P92 
P93 



pv-g + I^<yiy,a [cEy - AnByiT)] dv 



(14) 



E 



H 



\v\ 



\v\ 



+ 



p 



7 - 1 p 
7 P 



H 



2 7- 

(^'1, V2, V3) 



1 P 



The symbol, v = {vi, V2, fs), denotes the gas velocity and the symbol, g = {gi, g2, gs) 
the gravity. The gas is assumed to be ideal one of which the specific heat ratio is 7. 

2.3. Numerical Scheme 

Equations (2) and (5) can be expressed as 



(15) 

(16) 

does 



dt 



( E^ \ 

Fu,2 



c^P 
c^P 



{AtiB, - cE,) 

-C ((T^,a + <yu,s) Fy^i 
-C (cr^,a + <yv,s) Fy^2 
-C (0-^,a + <yu,s) F^^3 



(17) 



(18) 



Equation (18) has the same structure as that of the hydrodynamical equations in the conserva- 
tion form. Thus the Godunov-type method, which is often used for solving the hydrodynamical 
equations (see, e.g., Toro 2009), can be applied to Equation (18). 

In the Godunov-type method, the time evolution is evaluated based on the character- 
istics, i.e., the propagation speeds of signal. The characteristics of Equation (18) is rather 
complex and the computation of them takes much time. Gonzalez et al. (2007) obtained them 
by interpolating the table prepared rather than by computing them at each time step. 

We can avoid detailed computation of the characteristics by employing HLLE scheme, in 
which we need only the upper and lower limits on the characteristics (cf. Toro 2009). However, 
HLLE scheme gives us a too much diffusive solution if the upper and lower limits are taken to 
be the speed of light, ±c (Gonzalez et al. 2007). 

Ideal characteristics should be evaluated from the radiation fields in the adjacent cells 
across the surface on which the numerical fiux is evaluated. Figure 1 of Gonzalez et al. (2007) 
shows the characteristics for a given radiation field. The ideal characteristics should be appro- 
priately averaged ones when the radiation fields differ appreciably in the two adjacent cells. Roe 
(1981) obtained such average characteristics for the hydrodynamical equations. Similar average 
characteristics have not been obtained for the moment equation of the radiative transfer. 
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2.4- Kinetic Reconstruction of the Intensity 

We obtain the numerical flux not by computing the characteristics but by reconstructing 
the intensity consistent with the moments of the radiation field. When the intensity is expressed 
as 

/3. = /•^'^ (20) 

the moments as well as the closure relation, Equation (8), are consistent with the assumed 
radiation field. Although the closure relation cannot specify the intensity uniquely, Equation 
(19) has a distinctive feature: the entropy is minimum (Dubroca & Feugueas 1999). We 
can obtain the angular distribution by the Lorentz transform of an isotropic radiation field 
(Levermore 1984). 

We evaluate the flux, F^^i, across a cell boundary between two adjacent cells by using 
the reconstructed intensity. We call the two cells, left (L) and right (R), for later convenience. 
The cell surface is assumed to be normal to the i-th direction, i.e., {x^ — ccl) x Cj = 0. Then 
the intensity on the cell surface is evaluated to be 

W = . (21) 

[ h,R{n) [n-ei < 0) 

where Iu^i^{n) and /^^(n) denote the radiation fleld reconstructed in L and R cells, respectively. 
Equation (21) is based on the fact that photons transmit the surface from L to R when n Si > 
and from R to L otherwise. In other words, the intensity on the upwind side is identifled as 
that on the surface. Thus the energy flux, 

F;, = J ie,-n)i:{n)dn, (22) 

is expected to inherit the "upwind" nature and accordingly to be an alternative to the Godunov- 
type numerical flux. Similarly the momentum flux is evaluated to be 

P:,j = IJ {e,-n){e,-n)i:{n)dn, (23) 

The numerical flux, F*^, is expressed by an explicit function of E^i^, f^^i^, E^^^, and 
/^j^. For later convenience the numerical flux is decomposed into two components, 

^;.=4tL + ^ii' (24) 

i"i,U=/ (e,-n)/,,L(n)dn, (25) 
jnei>o 



.(-) 

The former is evaluated to be 



F^k= {e,-n)h,n{n)dn. (26) 

Jnei<o 
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7r/2 r2TT 



JO 



V^) COS 6* sin 6* dip 



^qu,h + 6/3^L COS^ ?/^^,L g,. ,L - /3^,L COS^^;^,Lg,.,L + 8/3,.,L COS ^^i^.L 



4(/3.\ + 3) 



\l/2 



9i^,L = (l - /^^Lsin^^Ai.,!) 
where the angular variables, 6', 9?, and ■?/', are defined to satisty 
(cj ■ n) = cos^ , 
(cj ■ n) = sin 6' cos ip , 
(cfc ■ n) = sin 6' sin ip , 
(Ci -/^l) = /3,.,L cos ?/;i.,L- 



(27) 

c^.,L, (28) 
(29) 

(30) 
(31) 
(32) 
(33) 



The symbols, and ek denote the unit vectors perpendicular to e^. The suffix, L, indicates 
that the variables are evaluated in cell L. We used the computer software, Mathematica, to 
obtain the above integral. Similarly we obtain the other half. 



TT /•27r 



7r/2 JO 



Iu,r{0, 'p) cos9 sm9d9 dip 



3g^,R + 6/3^ RCOs^?/;^,Rg^ R - Z^^^r cos^ z/;^,Rg^ - 8/3^,rcos?/^^,r 

4(/3^,R + 3) , 

Note that F^^i coincides with F^J}^ + F^ jj^ since ^ = 4/3^ (3 + /3^)-^ 



(34) 

cF,,R. (35) 




Fig. 1. Numerical flux obtained by the reconstructed intensity is shown as a function of / for 
(E, F = (1/c, 0, 0, /). The black solid and dashed curves denote F^^^ and Fi^\ respectively. The grey 
sohd and dashed curves denote Pzt'' and Pxt\ respectively. 

Figure 1 shows the numerical fluxes evaluated from the reconstructed intensity for 
(£", Fj., Fy, Fz) = (1/c, 0, 0, /). The black solid and dashed curves denote Fj"*") and 
Fj"*"^ as a function of /, respectively. Both Fj"*"^ and Fj"*"^ approaches to 1/4 in the limit of 
/ = (isotropic). They have the asymptotic forms of Fj+^ — )■ / and Fj+^ — y^(l — /)/8 in 
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the limit of / — ?■ 1. 

Similarly the radiation pressure is evaluated to be 



p* 

p(+) ^ 



p{-) 



(+) 



u,ii,h ' 
7r/2 p2-K 



pi-) 



l3^j^cos^ijjq~l + 3/3^,LCOs?/^j,,Lgi.,L + 4/32lCos2?/^^,l + 1 - 







2(/3^,L + 3) 
Iui{{9, cos^ 9 sm9d9 dip 

7r/2J0 

-/3^RCOs3?/^g"R - 3/3^.3 cos V',.,Rgj.,R + 4/3^ cosVi.,R + 1 - /^^r 



2(/3.',R + 3) 



P. 



(+) 



,(-) 



7r/2 ^27r 

h,Li9, ip) COS 6' sin^6' cos Lpd9 dip 











i/,ij,R 



TT /'27r 

7r/2 JO 



(36) 
(37) 

(38) 

(39) 

(40) 

(41) 
(42) 
(43) 
(44) 

(45) 



Iu,r{9, ip) cos9 sm'^9 cosipd9 dip 

= /3i.,j,R-^i,i,R- 

The values of Pi+) and Pi+) are denoted by the grey solid and dashed curves, respectively, in 
Figure 1. 

Remember that our numerical flux is obtained by integrating the intensity over the left 
and light hemispheres separately. Dubroca et al. (2003) proposed a similar idea named half 
space moment approximation. They integrated the radiative transfer equation over the half 
hemisphere and obtained two unknowns and two equations for each moment. Their equations 
are closed by two closure relations each of which is applied for each half hemisphere. Thus 
their scheme is different from ours in many respects, although the idea, integration over the 
half hemisphere is common. 

2.5. Inclusion of Absorption and Emission within a Cell 

In the previous subsection we implicitly assumed that the radiation field is uniform 
within a cell. This is not a good approximation when the cell under consideration itself is 
optically thick. The intensity differ appreciably on the cell surface from that at the cell center. 
Taking account of the absorption, emission, and scattering within a numerical cell, we modify 
Equation (24) into 



C^.L-^ilL + (1 - Cu,i,L] 



cE, 



4 



Cj.R-^itk ~ (1 ~ Cv,i,K, 



cE, 



(46) 
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'nu,i,L 


= exp 




— 




= exp 




= exp 



2 



(47) 
(48) 
(49) 
(50) 

Equation (46) denotes an approximation to the formal soultion of the radiative transfer in the 
hmit of Ax — )■ 0. On the other hand, in the hmit of of Ax — > oo, it describes the state in which 
the radiation is in the thermal equilibrium in each cell. 

Similarly Equations (36) and (41) are modified into 
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C^.R-PLi'Ir + (1 ~ C«,R; 



u,R 



6 



27r 

6C 



By^ij — 'nv,i,ljC,y,i,hPl,ij,l^ + Vv,i,RC,v;i,RPuAj,R 

We prove that our modified numerical flux is asymptotic preserving and reproduces the 
diffusion approximation in the limit of ftapAx — oo. For later convenience, we define the 
symbol 



(51) 
(52) 



AF,* 



E 



Axi 



Ax, 
X H — e, 



Ax,; 



(53) 



to denote the "divergence" of the numerical flux evaluated on the cell. Similary we define 



AP*.. 



V Axfc 



Axfc 



P, 



Axfc 

X — Cfc 



(54) 



to denote the divergence of the pressure tensor. The discretized Ml equations reduce to 

AF; = -K,p {cE, - AtxB,) , (55) 



(56) 



in radiative equilibrium. 

When Ax is large and the temperature difference is small. Equations (46), (51), and 
(52) reduce to 

F:^^^'n[B,{T^)-B,{T^)] (57) 



TT ■ 



dT 



Tt 



R 



27T 



Y[i?.(^L) + i?.(TR)] 



(58) 

(59) 
(60) 



respectively. 



10 



Substituting Equations (58) through (60) into Equations (55) and (56), we obtain 
c 

Air 1 dB^ dT 



Air 

E,^—B, (61) 

c 



''''' 3p Ky^a + K'V,s dT dXj ' 

where only the most dominant terms are taken into account for simphcity. These equations 
are equivalent to the diffusion approximation (see, e.g. Castor 2004, for frequency dependent 
diffusion approximation). Thus our numerical flux gives a good approximation in the optically 
thick limit. We can expect that our numerical flux gives a reasonable approximation at any 
optical depth. 

2. 6. Second Order Accuracy in Space 

The numerical flux given in the previous subsection is of the first order accuracy in 
space. A numerical flux of the second order accuracy in space can be obtained by applying 
Monotone Upwind-centered Scheme for Conservation Laws (MUSCL, see e.g., Hirsch 1990). 
MUSCL evaluates the physical states on the cell boundaries from the left and right hand sides 
by extrapolation. MUSCL applies a limiter such as minmod function in order to avoid spurious 
extrema from the values obtained by simple extrapolation. We need some additional cares to 
avoid unphysical extrapolation when applying MUSCL to Ml equations. Otherwise, the energy 
density can be smaller than the radiative energy flux divided by the speed of light {E < \F/c\). 
We use the formulae given in Appendix to achieve the second order accuracy in space. 

The numerical flux given in Appendix is constructed on the assumption that both emis- 
sion and scattering change linearly along the line from the cell center to the boundary. The 
emission and scattering on the cell boundary are evaluated by linear extrapolation with the 
minmod limiter, i.e., by MUSCL. Therefore the contribution of the emission and scattering 
within the cell is expressed as the linear combination of those evaluated at the cell center and 
boundary. Further details are given in Appendix. 

2.7. Time Evolution 

Using the notation given in the previous subsection, we can rewire the Ml equations in 
the form, 

= K,,aP {cE, - 47rfi,) - AF, , (63) 

dF 

= - {f^u,a + f^u,s) pcFyj - C^APyj . (64) 

We use two different methods to integrate Equations (63) and (64) in the first order accuracy 
in time depending on the problem. We use the forward difference of the first order in time in 
flash test in which we solve propagation of radiation from a sphere into vacuum. In the rest of 
test problems, we use a formal solution for forwarding the energy density and flux in time. 
The formal solution of of Equations (63) and (64) expressed as 
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V 



) 



(65) 



(66) 




where all the physical variables are evaluated at time, t. Thus our scheme is explicit in the 
sense and F^^j are obtained without iteration. Nevertheless, the formal solution allows to 
use a much longer time step than the simple forward difference in time when K^ apAx^ ^ 1. 
The time step should be smaller than At < Xji^Ky ^P) in the simple forward difference. When 
Ky^aP^x ^ 1, the constraint is serious since it is much shorter than the time for propagation 
of signal, As/c. 

We integrate the hydrodynamical equations by an explicit scheme. The time step is 
taken to be the minimum of the thermal and hydrodynamical timescales. 

A solution of the second order accuracy in time can be obtained by a two step Runge- 
Kutta method. We used the average of time derivatives evaluated at t = to ^"^^ + when 
integrating equations from t = to to to + At. 

3. Monochromatic Test Problems 

3.1. Shadow Test 

We performed a shadow test to illustrate the effects of including absorption within 
the cell. We exposed uniform monochromatic radiation to a square absorber of np = 50. 
Scattering and emission are neglected for simplicity. The computation box covers the square 
of < x < 12 and < y < 6. The absorber occupies the square region of 2 < x < 3 and 
< y < 2. The spatial resolution is Ax = Ay = 0.1. We imposed the uniform radiation, 
{E, Fx, Fy) = (1.0, 0.999 c, 0), from x = 0.0. Reflection boundary is applied to the upper 
and lower boundaries of y = 6 and 0. The outgoing boundary is placed on the right boundary, 
X = 12.0, so that no radiation enters from the boundary. We realized the outgoing boundary 
by vanishing the flux from outside the boundary. It is a virtue of the kinetic flux that the 
outgoing boundary is easily constructed. We obtained the equilibrium state by solving the 
time dependent Ml model with the time step. At = 0.5 Ax/c. 

Figure 2 shows the result of the shadow test in the equilibrium state. The numerical flux 
of the first order accuracy in space is used in the upper panel while that of the second order 
accuracy is used in the lower panel. The brightness denotes the energy density in logarithmic 
scale. Note that the energy density drops very sharply behind the left side of the absorber in 
the time step. The arrows denote F/{cEy). 

Very weak radiation shines in from the upper right corner behind the absorber because 
the incident radiation is not perfectly beamed {F/{cE) = 0.999). FWHM of the beam is 
^fwhm/2 = l.°58, since it is evaluated to be 6'fwhm/2 — 0.870-^/1 — / from Equations (19) 
and (20) for 1 - /3 ~ 2(1 - /) < 1. The intensity decreases by a factor 100 at ^ = 5.°33. The 
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inclination of the contour of log E = — 2 is 6° and 5° from the horizontal line in the solutions 
of the first and second order accuracy, respectively and constant with the intensity distribution 
of / = 0.999. The solutions are different only in dark area of log E < — 2. The contour of 
log E = — 3 is more inclined in the solution of the first order accuracy than in that of the 
second order accuracy. 
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Fig. 2. 2 Dimensional Shadow Test. The upper and lower panels show results obtained by the numerical 
flux of the first order accuracy in space and that of the second order, respectively. The dashed line denotes 
the absorber. The brightness denotes the radiation energy density in the logarithmic scale. The solid 
curves denote the contours of log£^ = —3,2, and — 1. See text for more details. The arrows denote the 
vector, Fl\E\ 



3.2. Flash Test 

We performed the following test for studying the propagation of radiation in vacuum. 
The initial radiation field at t = is set to be 

(B, f^, f,, = I (1. c/. 0. 0) ,fH<r„ 

^ " ^ I (0, 0, 0, 0) if |r| > ro ^ ^ 

so that uniform radiation is confined within the sphere of r < tq. No absorption and emission 
are assumed in this test. Then the radiation is expected to be confined in the spherical shell of 
Tq — ct < r < ro + ct. 

Figure 3 denotes the radiation energy density for / = 0.0 and ro = 0.5 at t = 6.0/c. 
The spatial resolution and time step are taken to be Ax = 0.1 and At = 0.3 Ax/c, respectively. 
The upper panel denotes the result obtained by the simple HLLE flux, while the lower panel 
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does that obtained with our numerical flux. Both of them denote the solutions of the first order 
accuracy in space and in time. The brightness and the contours denote logE. It should be 
noted that HLLE does not work at At = 0.35 Ax/c, while our method does at At = 0.5 Ax/c. 
A longer time step can be taken safely if we apply the kinetically reconstructed numerical flux. 




X 



Fig. 3. The result of flash test for / = 0. The upper panel shows the solution by HLLE scheme and 
the lower panel does that by our numerical flux. Both the panels denote the solutions of the first order 
accuracy in space and in time. 

If we solve the flash problem perfectly by taking the full angular dependence, the radi- 
ation should be confined an expanding spherical shell of ct — tq < r < ct + tq. It should be 
isotropic around the origin and the flux should be radial. 
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Our model shows a higher contrast than the simple HLLE model. However a clear 
anisotropy of numerical origin appears in our model. This anisotropy is due to the fact that 
the change in the flux direction in a cell is not taken into account. 

Figure 4 is the same as Figure 3 but for the solutions of the second order accuracy 
in space and in time. The anisotropy has been removed in the solutions of the second order 
accuracy. The low contrast of the simple HLLE model is also improved, although the contrast 
is still higher in our model. 



logE 




X 



Fig. 4. The same as Fig. 3 but for the solutions of the seeond order accuracy in space and in time 

Figure 5 is the same as Figure 3 but for / = 0.9. When / is close to 1.0, the forward 
radiation is much stronger than the backward one. However, the shell of the high radiation 
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energy density should expand spherically as in case of / = 0, if we take the full angular 
dependence of the radiation. It expands as expected in the lower panel (kinetically reconstructed 
numerical flux), while the center of the sphere shows a spurious shift in the upper panel (HLLE). 
This is due to the characteristic speeds of waves in Ml model equations. Their absolute values 
are smaller than the speed of light and all the characteristic speeds have the same sign when 
/ > 0.69 (see e.g., Figure 1 of Gonzalez et al. 2007 for the characteristic speeds as a function 
of /). Thus the center of radiation shell shifts rightward in the HLLE solution. 




Fig. 5. The same as Fig. 3 but for / = 0.9. 

One might think that the shift would be the nature of Ml model equations and should 
be reproduced in numerical simulations. However, we should realize that the characteristic 
speed changes according to that in / = F/{cE). The ratio, /, decreases on the left edge of 
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the radiation sphere, since radiation moving rightward has a higher / than the initial value. 
Once it decreases to the critical value on the edge, the radiation begins to propagates leftward 
and it approaches to / = — 1. It should take only an instant for the left edge of the radiation 
sphere to propagate at the speed of — c, if our spatial resolution is extremely high. However, it 
takes a few time steps for / to decrease to the critical value in a cell close to the left edge. The 
propagation to rightward is delayed a few time steps in the HLLE model. 

Figure 6 shows the second order accurate solutions for the flash test of / = 0.9. The 
bright shell is more sharply captured in both the second order accurate solutions. The contrast 
is still higher in our model than in the HLLE model. The central dark hole is shifted leftward 
in the HLLE model. 
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Fig. 6. The same as Fig. 5 but for the solutions of second order accuracy in space and in time. 
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4. Application to Irradiated Protoplanetary Disks 



Young stars are often associated with gaseous disks called protoplanetary disks (see, e.g, 
a review by Williams & Cieza 2011). They are irradiated by the radiation from the central stars 
to shine at various wavelengths. They reflect optical and near infrared stellar lights and emit 
mid and far infrared ones. Both absorption and scattering are dominated by dusts, although the 
optical properties remain somewhat unknown. It is essential to take account of the frequency 
dependent opacity when modeling protoplanetary disks. Optical and near infrared radiation 
heat up the disks from outside while mid and far infrared emission cool down them from inside. 
The mid and far infrared flux from inside balances the optical and near infrared one from 
outside in equilibrium. 

We apply our numerical method to an irradiated protoplanetary disk as a test for multi- 
color problem. We use the opacity table of Draine (2003) in which Ki^q, k^^s, and (cos 6') are 
given as a function of the wavelength. Applying the sphne fit to the table, we obtained the 
values at the wavelengths, 

loglY^^) = 0.02 m, (68) 
y 1.0 /im J 

where m is a integer in the range of —50 < m < 150. Figure 7 shows K,,y^a, i^u,s, and 
i^u,a + i^u,s{^ — {cos 6)) as a function of A. The opacity is obtained under the assumption that 
the dust occupies 1% of the total mass. In other words, we did not take account of sedimentation 
of dust for simplicity. 




-^ i j 

-10 12 3 

log X 



Fig. 7. The opacity used in our application to irradiated protoplanetary disks. The curves denote, k, 
Ki^^s, and K^^a + 'Sj/.s (1 ^ (cos(?)) as a function of the wavelength, A = c/v. 



We use two sets of Ml model; one denotes the direct radiation from the central star 
and the other does scattered by and reemitted from the protoplanetary disks. Then radiation 
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energy density and flux are expressed as 

E, = El + E'l, (69) 

F, = n + F'J, (70) 

where the symbols with prime denote the values of the direct stellar radiation and those with 
double prime do the total of scattered radiation and emission from the disk. We apply the 
closure relation separately for the two components. The Ml model equations are expressed as 
f)F' ^ dF' 

^ + E ^ = - + K.,s (1 - (cos^)] pcEl . (71) 

i=l '* 

^ + E ^ = (1 - (cos^))] pcEl - K.,,p [cE: - 47ri?.(T)] , (72) 
OF' 3 

^ + En., = - K,a + n^,s (1 - (COS^)] , (73) 

dF" 3 

+ c'E^:.. = - + ^u,s{i - {oose))]cFr. (74) 

^ i=i 

The separation of the direct stellar light from the rest radiation avoids spurious beam collision 
on the disk surface. The radiation energy densities are evaluated at the wavelengths given by 
Equation (68), i.e., 201 bands in the range of 0.1 yum < A < 1 mm. Thus our model has the 
spectral resolution of AA = 4.61 x 10~^ A. This spectral resolution is good enough to study 
the spectral features of dust opacity. 

We solved the above moment equation of radiation and the hydrodynamical equations 
simultaneously. In this paper we restrict ourselves to the disk in equilibrium in which the 
emission from the disk balances the heating by irradiation. We ignored self gravity of the disk 
and viscous heating by accretion for simplicity. The radial component of the gravity is assumed 
to be balanced with the centrifugal force due to the disk rotation. 

We assume that the central star has the mass and radius, M^, and i?*, respectively. The 
stellar radiation is assumed to be the black body of Teg. We made both one and two dimensional 
models of the protoplanetary disk with the numerical flux of the first order accuracy in space. 
They are described in subsections 4.1 and 4.2, respectively. 

4-1. ID Model Based on the Grazing Recipe 

Our ID model describes the vertical structure of the protoplanetary disk at a given 
radius, R, from the central star. The stellar radiation, E[„ is evaluated to be 



TT fR. 



EUr, z) = -^-^) i?,(Teff)exp 



2 



T,AZ] 



a 



(75) 



T,AZ 



POO 

[«v,a + i^y,s (1 - (cos6'))] / p(r, z') dz' , (76) 

J z 
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according to the grazing angle recipe (Chiang & Goldreich 1997). Here if* denotes the height of 
the 'photosphere' at which the stellar radiation is attenuated by a factor of e~^. We evaluated 
the photosphere at A = 0.302 fim. We evaluated d{H^/ R)/dR consistently by solving the 
vertical structure at a slightly different radii, 0.891 R and 1.122 R. 

We obtained the steady state solution by integrating the Ml model equations and equa- 
tions for hydrostatic balance simultaneously. We used the Lagrangian coordinate in this ID 
model. 

Figure 8 shows the model spectra at i? = 50, 100, and 200 AU from the star of 
M^, = 2.0 Mq, -R* = 2.5 -Rq, and T^s = 9,500 K. These parameters are taken to be similar 
to those of AB Aur (van den Ancker et al. 1997) for which inner hole and spirals structure are 
seen in optical (Grady et al. 1997) and near infrared (Fukagawa et al. 2004; Hashimoto et al. 
2011). 

The surface density is assumed to be 

The model spectra are consistent with those obtained by D'Alessio et al. (1998); DuUemond, 
Dominki, & Natta (2001), who solved the angle dependent radiative transfer with the Monte 
Carlo simulation. The ratio of the grazing angle to the disk aspect is d/dhiR\ln{H^,/ R)] = 0.251 
at R = 100 AU, which is close to the standard value, 2/7 = 0.282 (Chiang & Goldreich 1997). 
The value depends rather on the opacity and surface density distribution. The standard value 
was obtained by a simplified analytical model. The difference is not numerical. 




log X ( nm) 

Fig. 8. Model spectra for a protoplanetary disk irradiated by the star of M = 2.0 Mq, R ~ 2.5 i?*, 
and Teff = 9500K. Each curve denotes the flux from the disk surface at the designated radius. 



Figure 9 confirms that the heating is balanced with the cooling correctly in our solution. 
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Each curve denotes the energy flux, vFl^ and z/F^'^ a given height as a function of the 
wavelength. The value is taken to be positive when the energy flows outward from the disk. The 
solid curves denote F^^ (negative) and F^^^ (positive) at the surface (not the photosphere). The 
dashed curve denotes F"^ at the level above which the disk has the surface density, 3.9 x 10~^ 
g cm~^. The dash dotted curve does that at S = 3.5 x 10~^ g cm~^. It is clearly shown that 
the net flux vanishes at any hight. The mid infrared is the main heating source in the layer 
below S > 3.5 x 10"^ gm cm~^. 



+ out 



E 
o 

D3 
> 



in 




-600 



surface 

E = 3.9x1 " g cm"^ 
E = 3.5x1 0"^ g cm"^ 



1 

log l ( \im) 



Fig. 9. The wavelength dependent energy flux. The surface density denotes that above the layer. 



4-2. 2D Axisymmetric Model 

ID model assumes implicitly that the surface density changes gradually with the radius. 
However some protoplanetary disks may have holes and the surface density may rise sharply at 
some radius. The transition disks are thought to be the case (see, e.g., a review by Williams & 
Cieza 2011 and the references therein). We construct a two dimensional model of a transitional 
disk assuming the symmetry around the the axis and that on the mid plane. 

Following Honda et al. (2012), who made spectral energy distribution (SED) model for 
a transitionl disk around HD169142, we made a model in which the surface density is expressed 
as 



InE = - <{ In(SinSout) + In 



-'out 



tanh 



rin 



m 



(79) 



where Sin, ^out, E, and tq are model parameters to be chosen. We obtained the radiation and 
gas in the region of 40 AU < r < 160 AU and 1^1 < 80 AU, with the resolution of Ar = 0.3 AU 
and Az = 0.4 AU in the cylindrical coordinates. We placed the reflection boundary on Z = 
and the outgoing boundary on z = 80 AU. The incoming flux from R = 40 AU was fixed. 
The flux incoming from r = 160 AU was assumed to balance with the outgoing one in the disk 
and to vanish outside the disk. The mass, radius, and effective temperature of the central star 
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are assumed to be M = 2.0 Mq, = 1.25 x 10"^ AU, and T^s = 9,000 K, respectively. 

Figure 10 shows the density and temperature distribution in equihbrium for ro = 100 AU, 
r = cx). Sin = 3.5 X 10-3 g cm-3, and Sout = 3.5 x lO^^ g cm-^. The wall at r = 100 AU is 
heated up to T ^ 140 K and hence expands more. The grayness denotes the temperature and 
the curves denote the contours of logp in the interval of Alogp = 0.5. Note that we solved the 
vertical hydrodynamic balance consistently while the vertical density distribution is given in 
Honda et al. (2012). The radiation from the wall heats up the inner disk by 15 K in the range 
of 70 AU< r < 100 AU compared with the model without the wall, i.e., S = Sin {r/ro)~^. 




Fig. 10. The temperature and density distribution in the model having the sharp rise in the surface 
density at r = 100 AU. The greenness denotes the temperature and the curves denote the contours of 
logp in the interval of Alogp = 0.5. 

Figure 11 shows the total energy density, E^, at A = 0.316 /im, 1.58 yum, 20 fim, and 
501 yum from top to bottom. The color denotes logEi, in unit of erg cm"'^ Hz"^ as indicated in 
the left bar. The arrows denote the vector F^/E^. 

At A = 0.316 /im, stellar radiation absorbed or scattered on an upper layer of low 
density and does not penetrate into protoplanetary disks. Near infrared radiation penetrate a 
little deeper interior of the disk but does not reach the mid plane. The upper part of the wall 
is bright at A = 20 /im. The wall is heated up by stellar light and emits mid and far infrared 
radiation. The disk is more transparent at a longer wavelength as shown in the bottom panel 
of Figure 11. The direction of energy flux also depends on the wavelength. The energy flows 
from the central star in the optical and near infrared, while it flows from the wall and upper 
surface in the mid infrared and from the disk in the far infrared. 

We obtained simulated images of the protoplanetary disk by integrating Equation (1) 
along the line of sight. The source terms are evaluated from T and E^, obtained by our Ml 
model. The upper panels of Figure 12 show the simulated images at A = 1.58 /xm (left, 
H-band) and 20 /im and (right, Q-band). The line of sight is assumed to inclined by 15° from 
the axis normal to the disk. Both the images show bright rings at r ~ 100 AU. These images 
are similar to those of Honda et al. (2012) who solved the radiative transfer by the Monte Carlo 
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simulation. 

We made another model by assuming F = 20 while keeping the other model parameters 
unchanged. The temperature and density distributions are shown in Figure 13. The density 
change around the wall is smooth and more likely than that shown in Figure 10. 

Figure 14 is the same as Figure 11 but for F = 20. The result is almost the same but 
the wall boundary is less sharp as expected. The photospheres are located just behind the wall 
in the model of F = oo in a broad range of the wavelength. However, the location of the 
photosphere depends on the wavelength in the model of F = 20. 

The lower panels of Figure 12 is the same as the upper panels but for F = 20. The 
bright ring are broader in the model of F = 20 than in that of F = oo. The peak brightness 
is lower in the model of F = 20. This is because the wall is appreciably inclined when F = 20 
(see Figure 13). The ring corresponds the region in which the surface density increases with 
the radius. In other words, the wall is seen as the bright ring in the image. Our model will 
serve to evaluate the surface density distribution from observed images. 

5. Discussions 

As demonstrated in the previous sections. Ml model works well both in vacuum and in 
optically thick media if absorption is taken into account properly in the numerical flux. If we 
use the formal solution for computing time evolution, the time step can be as large as Ax/ (2c) 
irrespectively of the opitacal depth. Thanks to this improvement, we succeeded in applying Ml 
model to the protoplanetary disks. They are optically thick in optical bands while optically in 
far infrared bands. They are heated by the optical and near infrared stellar lights and cooled 
by mid and far infrared emission. Thus it is important to take account of both optically thin 
and thick radiation simultaneously. Our numerical scheme will expand the applicability of Ml 
model. 

Ml model can be used also for neutrino transfer. Neutrino is a major coolant in compact 
objects such as neutron stars and black holes. It should play an important role in dynamics 
of core collapse supernovae and in gamma ray bursts (see, e.g., a review by Mezzacappa 2005 
and the references therein). Both core collapse supernovae and gamma ray burst sources are 
thought to be highly anisotropic and their dynamics should be studied ideally based on three 
dimensional numerical simulations. It is also essential to take account of the neutrino energy 
in the simulations. Neutrino opacity is proportional to the square of the neutrino energy. Low 
energy neutrinos are easy to leak while high energy ones not. Thus Ml model is a reasonable 
choice for numerical simulations of these objects. It reduces the computation cost by lowering 
the angular resolution. Yet it can express a shadow and beamed radiation. 

It should be noted that Ml model can solve the propagation of a flash by an explicit 
manner. Radiation at the next time step depends only on those in the neighboring cells. We 
need neither ray tracing nor global iteration. In other words. Ml model does not require global 
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communication for proceeding a time step. This is beneficial for massive parallel computation, 
since global communication between processors are often a bottle neck. 

In compact objects, the gas has a high temperature and the sound speed is comparable 
to or only by a factor of ten smaller than the speed of light. Thus it is not serious that the 
time step is restricted to the light propagation time, 0.5 Ax/c; a similar small time step is 
required from the hydrodynamical simulations if we integrate them explicitly. If we solve both 
radiation and hydrodynamics explicitly, we can take account of heating by a flash of neutrino. 
The compact sources may have a highly variable luminosity. 

Ml model may be applied to dynamics of non-relativistic objects in which the sound 
speed is much slower than the speed of light. We can reduce the speed of light propagation 
in Ml model to prolong the time step. If we replace c by d in Equations (65) and (66), the 
propagation speed is reduced to c'. We expect that the reduction does not affect the result 
seriously from the following reason. Both hydrodynamical and thermal timescales are much 
longer than the time for light propagation in most non-relativistic objects. Thus the light 
propagation speed is assumed to be infinite in some radiation hydrodynamics. The results are 
valid since the light speed is much faster than other speeds of waves and we can neglect the 
difference between the real light speed and infinity. We think that we can neglect the difference 
between c and d as long as d is much larger than other speed of wave propagation. Our 
idea is based on the heuristic experiment by Hotta et al. (2012). They performed numerical 
simulations of convection in the Sun by reducing the sound speed artificially. The velocity of 
material convection is much lower than the sound speed. Thus, most of the simulations thus far 
apply the anelastic approximation in which the sound speed is infinitely large. However, Hotta 
et al. (2012) demonstrated that nearly the same results are obtained even when the sound speed 
is reduced artificially as far as the reduced sound speed is still much faster than the convection 
velocity. Their experiment suggests us that we can reduce the light speed artificially without 
loss of quality. 

In summary. Ml model has potential applicability to many problems including proto- 
planety disks, core collapse supernovae, and gamma ray bursts. 
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Appendix. Numerical Flux of the Second Order Accuracy in Space 



First we introduce a new variable, 



E' 



El-\ — 



(Al) 



for later convenience. We evaluate cell boundary by extrapolation with the 

minmod limiter and obtain 

1 



p/(L)* 



E' 



2 ^'.i+i/s' 



e: 



'(R)* 

i^,i+i/2 



p/ J- . p/(R) 
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!^J+l/2 



X 



(A2) 
(A3) 

(A4) 

(A5) 
(A6) 
(A7) 
(A8) 

where sgn denotes the sign function. Similarly we evaluate -Fa;,i/,j+i/2 on the cell boundary to 
obtain 
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(A12) 
(A13) 
(A14) 
(A15) 

We obtain F^}*^^i^, F^u*j+i/2^ ^iS+1/2 and i^iJj+1/2 by the same procedure. The energy 
density on the cell boundary is evaluated to be 
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Introduction of E'^ guarantees that the energy density is larger than the energy flux divided 
the speed of light on the cell boundary. 

The energy density, E^^^-^^^, can exceed {?>/2)Ey ,j. If we would apply MUSCL to E^, it 
can not so high. To avoid such a high energy density, we reduce both the energy density and 
flux by multiplying the factors. 
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Accordingly we obtain 
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We use E^j^-^^2 ^-^d -^"[^^-+1/2 for computing the energy flux from cell j to j + 1, and E^j_^_-^^^2 
and F)^jj^ii2 fo^' that from j + 1 to j. 

Emission and scattering within a numerical cell are included in the second order numer- 
ical flux by the following procedure. First we evaluate the the blackbody emission on the cell 
boundary, -B^^^_,_i/2 ^'^'^ -^1^+1/2 ^"7 the MUSCL approach. 
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We assume that Bi, is a linear function of the optical depth between the cell center and boundary. 
Then we obtain 



K, 



■,a.,jPjBu{x ) dx — -B^ j_|_^/2 
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1 - e-^'-J 



w 



(A28) 



where 



Wu,j = K^,aP^Xj. (A29) 

The symbol Axj denotes the cell width. Similarly we can evaluate the scattering within the 
cell. 

By taking the emission and scattering within the cell, the numerical flux of the second 
order accuracy is expressed as 
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11. Total energy distribution, Ey(r^ z) at 0.316 /im, 1.58 yum (H-band), and 20 yum (Q-band), and 
/im from top. 
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Fig. 12. Simulated images of the protoplanetary disks with a wah at H and Q bands. The upper panels 
are for the model with a sharp edge (F = oo) while the lower panels are for the model with a smooth 
transition (F = 20). The left panels denote the brightness in H band in the logarithmic scale while the 
right panels do those in Q band. 
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Fig. 14. The same as Figure 11 but for the model having F = 20. 
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